*Title: Delivering remote learning using a low-tech solution: Evidence from an RCT during school closures  
*Author: Liang Choon Wang, Michael Vlassopoulos, Asad Islam, Hashibul Hassan 

*These codes generate Figure B3-B8 of the Appendix

***************************************************************************************
*Stata preparation
***************************************************************************************
*ssc install schemepack, replace
*net install grc1leg2.pkg, replace
graph set window fontface "Times New Roman"
set scheme white_tableau

***************************************************************************************
*Figure B3. Distribution of total number of lessons completed by caregiver-child dyads, by module
***************************************************************************************
cd ..
pwd
use "Data\IVR_Data.dta", clear
cd "Output"

preserve
* keep only treated
keep if any_treat==1

* generate three modules per person
expand 3
bysort CHILD_ID: gen Module = _n
lab def lblmodule 1 "Literacy" 2 "Numeracy" 3 "Leadership", replace
lab val Module lblmodule

gen mod = "lit" if Module==1
replace mod = "num" if Module==2
replace mod = "lead" if Module==3

* drop leadership module for standard group
drop if Module==3 & treatment_arms==1

* expand to 32 different lesson numbers (from 0, access, 1, ..., 30)
gen LessonNumber = .
foreach v in lit num lead {
	replace LessonNumber = totlessons_`v' if mod=="`v'"
	replace LessonNumber = 0 if evercomp_`v'==0 & mod=="`v'"
	replace LessonNumber = 31 if totlessons_`v'==. & evercomp_`v'==1 & mod=="`v'"
}
drop if evercomp_lit==.
gen Num = 1
* collapse by Module and LessonNumber for graphing
collapse (count) Num, by(Module LessonNumber)
lab def lblnum 0 "0" 1 "1" 2 "2" 3 "3" 4 "4" 5 "5" 6 "6" ///
				7 "7" 8 "8" 9 "9" 10 "10" 11 "11" 12 "12" 13 "13" ///
				14 "14" 15 "15" 16 "16" 17 "17" 18 "18" 19 "19" ///
				20 "20" 21 "21" 22 "22" 23 "23" 24 "24" 25 "25" ///
				26 "26" 27 "27" 28 "28" 29 "29" 30 "30" 31 "A"
lab val LessonNumber lblnum

graph bar Num, over(Module) over(LessonNumber, label(labsize(vsmall))) ///
asyvars ytitle(Number of caregiver-child dyads, size(small)) ///
subtitle(Number of lessons completed, size(small) position(6)) ///
ylabel(0(100)500, ang(h) labsize(small)) ///
bar(1, fcolor(red*.5) lcolor(black)) ///
bar(2, fcolor(blue*1) lcolor(black)) ///
bar(3, fcolor(yellow*1) lcolor(black)) /// 
legend(rows(1) ring(0) symxsize(8) position(12) size(small) region(lstyle(none) fcolor(none))) ///
graphregion(color(white)) name(FigB3, replace)

restore

graph export "Figure B3.emf", replace

***************************************************************************************
*Figure B4. Completion by caregiver-child dyads, by module and lesson number
***************************************************************************************
cd ..
pwd
use "Data\IVR_Data.dta", clear
cd "Output"

preserve
* keep only treated
keep if any_treat==1

* generate three modules per person
expand 3
bysort CHILD_ID: gen Module = _n
lab def lblmodule 1 "Literacy" 2 "Numeracy" 3 "Leadership", replace
lab val Module lblmodule

gen mod = "lit" if Module==1
replace mod = "num" if Module==2
replace mod = "lead" if Module==3

* drop leadership module for standard group
drop if Module==3 & treatment_arms==1

* expand to 30 different lesson numbers (1, ..., 30)
expand 30
bysort CHILD_ID Module: gen LessonNumber = _n

* drop lesson 17+ for leadership
drop if LessonNumber>15 & Module==3

* generate indicator for lessons listened in each module
gen NumberofListener = .

* fix reported 1 to 30 to 2 to 31 for figure labeling
forvalues i = 1(1)30 {
	replace NumberofListener = 1 if complete_lit_`i'==1 & Module==1 & LessonNumber==`i'
	replace NumberofListener = 1 if complete_num_`i'==1 & Module==2 & LessonNumber==`i'
}
forvalues i = 1(1)15 {
	replace NumberofListener = 1 if complete_lead_`i'==1 & Module==3 & LessonNumber==`i'
}

* collapse by Module and LessonNumber for graphing
collapse (sum) NumberofListener, by(Module LessonNumber)
drop if (LessonNumber==0)
lab def lblnum 1 "1" 2 "2" 3 "3" 4 "4" 5 "5" 6 "6" ///
				7 "7" 8 "8" 9 "9" 10 "10" 11 "11" 12 "12" 13 "13" ///
				14 "14" 15 "15" 16 "16" 17 "17" 18 "18" 19 "19" ///
				20 "20" 21 "21" 22 "22" 23 "23" 24 "24" 25 "25" ///
				26 "26" 27 "27" 28 "28" 29 "29" 30 "30"
lab val LessonNumber lblnum

graph bar NumberofListener, over(Module) over(LessonNumber, label(labsize(vsmall))) ///
asyvars ytitle(Number of caregiver-child dyads, size(small)) ///
subtitle(Lesson number, size(small) position(6)) ///
ylabel(0(200)1000, ang(h) labsize(small)) ///
bar(1, fcolor(red*.5) lcolor(black)) ///
bar(2, fcolor(blue*1) lcolor(black)) ///
bar(3, fcolor(yellow*1) lcolor(black)) /// 
legend(rows(1) ring(0) symxsize(8) position(2) size(small) region(lstyle(none) fcolor(none))) ///
graphregion(color(white)) name(FigB4, replace)

restore

graph export "Figure B4.emf", replace

***************************************************************************************
*Figure B5. Total hours of audio lessons accessed, by program week
***************************************************************************************
cd ..
pwd
use "Data\IVR_Data_Usage.dta", clear
cd "Output"

twoway (line tothours week, xline(11) ///
		ylabel(0(100)800, ang(h)) ///
		ytitle("Total hours of lessons accessed") ///
		xlabel(1(1)15) xtitle("Program week")), ///
		name(FigB5, replace)

graph export "Figure B5.emf", replace

***************************************************************************************
*Figure B6. Assessment-test score distribution, by treatment 
***************************************************************************************
cd ..
pwd
use "Data\IVR_Data.dta", clear
cd "Output"

*Histogram
summarize endline_total_score if treatment_arms == 3
local m1 = r(mean)
local s1 = r(sd)

summarize endline_total_score if treatment_arms != 3
local m2 = r(mean)
local s2 = r(sd)

twoway 	(histogram endline_total_score if treatment_arms==3, start(0) width(3)) ///
					(histogram endline_total_score if treatment_arms!=3, start(0) width(3) fcolor(none) lwidth(.2) ) ///
					(function normalden(x,`m1',`s1'), range(endline_total_score)) ///
					(function normalden(x,`m2',`s2'), range(endline_total_score) lpattern(dash)), ///
					legend(order(1 "Control" 2 "Treatment" 3 "Normal Density (Control)"  4 "Normal Density (Treatment)") row(2) position(11) ring(0)) /// 
					xtitle(Total score) ytitle("Density, by treatment group") name(histogram, replace) title({bf:Panel A: Histogram}, size(medsmall))

*Percentile plot
drop if endline_completed==0

twoway__histogram_gen endline_total_score if treatment_arms == 1, percent  gen(ht1 xt1)
twoway__histogram_gen endline_total_score if treatment_arms == 2, percent  gen(ht2 xt2)
twoway__histogram_gen endline_total_score if treatment_arms == 3, percent  gen(hc xc)

keep xt1 ht1 xt2 ht2 xc hc

drop if xc ==. 
drop if hc==. 

g cum_hc= sum(hc)
g cum_ht1= sum(ht1)
g cum_ht2= sum(ht2)

line cum_ht1 cum_ht2 cum_hc, name(percentile, replace) msymbol(none) connect(l) xtitle(Percentile of control group) ytitle(Percentile of treatment groups) ///
title({bf:Panel B: Percentile to percentile plot}, size(medsmall)) legend(order(1 "Standard" 2 "Extended")  position(10) rows(1) ring(0)) ylabel(0 (10) 100) xlabel(0 (10) 100) || line cum_hc cum_hc, color(gray)

graph combine histogram percentile, cols(1) name(Figure_B6, replace)
graph display, xsize(6) ysize(8)

graph export "Figure B6.emf", replace

***************************************************************************************
*Figure B7. Dosage responses in standardized literacy and numeracy scores
***************************************************************************************
cd ..
pwd
use "Data\IVR_Data.dta", clear
cd "Output"

****Literacy against total lessons completed ****
preserve
drop if evercomp_lit==0 & any_treat==1
drop if endline_completed~=1
xi: reg endline_literacy_score_std i.cbin_lit, cl(VILLAGE_ID)
gen semean = _se[_cons] if any_treat==0
forvalues i = 1(1)4 {
	replace semean = _se[_Icbin_lit_`i'] if any_treat==1 & cbin_lit==`i'
}
collapse (mean) mean = endline_literacy_score_std semean, by(cbin_lit treatment_arm)
gen ll = mean - semean*invttail(89, 0.025)
gen ul = mean + semean*invttail(89, 0.025)

replace cbin = cbin - 0.05 if treatment_arms==1
replace cbin = cbin + 0.05 if treatment_arms==2

twoway(scatter mean cbin if treatment_arms==1, msymbol(square)) ///
	(scatter mean cbin if treatment_arms==2, msymbol(square)) ///
	(scatter mean cbin if treatment_arms==3, msymbol(round) ///
	xlabel(0 "0" 1 "1-19" 2 "20-29" 3 "30" 4 "Some", labsize(small)) xtitle("Number of literacy lessons completed") ///
	ylabel(-0.2(0.2)1.2, ang(h) labsize(small)) ///
	ytitle("Literacy score (standardized)")) ///
	(rcap ll ul cbin if treatment_arms==1, blwid(thin) blcolor(black) msize(medium)) ///
	(rcap ll ul cbin if treatment_arms==2, blwid(thin) blcolor(black) msize(medium)) ///
	(rcap ll ul cbin if treatment_arms==3, blwid(thin) blcolor(black) msize(medium)), ///
	legend(order(1 "Standard" 2 "Extended" 3 "Control" 4 "95% CI") cols(4) colgap(2) position(6) size(small)) ///
	graphregion(color(white)) graphregion(color(white)) ///
	name(FigureB7a, replace)
restore


****numeracy against total lessons completed ****
preserve
drop if evercomp_num==0 & any_treat==1
drop if endline_completed~=1
xi: reg endline_numeracy_score_std i.cbin_num, cl(VILLAGE_ID)
gen semean = _se[_cons] if any_treat==0
forvalues i = 1(1)4 {
	replace semean = _se[_Icbin_num_`i'] if any_treat==1 & cbin_num==`i'
}
collapse (mean) mean = endline_numeracy_score_std semean, by(cbin_num treatment_arm)
gen ll = mean - semean*invttail(89, 0.025)
gen ul = mean + semean*invttail(89, 0.025)

replace cbin = cbin - 0.05 if treatment_arm==1
replace cbin = cbin + 0.05 if treatment_arm==2

twoway(scatter mean cbin if treatment_arm==1, msymbol(square)) ///
	(scatter mean cbin if treatment_arm==2, msymbol(square)) ///
	(scatter mean cbin if treatment_arm==3, msymbol(round) ///
	xlabel(0 "0" 1 "1-19" 2 "20-29" 3 "30" 4 "Some", labsize(small)) xtitle("Number of numeracy lessons completed") ///
	ylabel(-0.2(0.2)1.2, ang(h) labsize(small)) ///
	ytitle("Numeracy score (standardized)")) ///
	(rcap ll ul cbin if treatment_arm==1, blwid(thin) blcolor(black) msize(medium)) ///
	(rcap ll ul cbin if treatment_arm==2, blwid(thin) blcolor(black) msize(medium)) ///
	(rcap ll ul cbin if treatment_arm==3, blwid(thin) blcolor(black) msize(medium)), ///
	legend(order(1 "Standard" 2 "Extended" 3 "Control" 4 "95% CI") cols(4) colgap(2) position(6) size(small)) ///
	graphregion(color(white)) graphregion(color(white)) ///
	name(FigureB7b, replace)
restore

graph combine FigureB7a FigureB7b, ///
	iscale(.8) col(2) name(FigB7, replace) graphregion(color(white))

graph export "Figure B7.emf", replace

***************************************************************************************
*Figure B8. Differential Item Functioning (DIF). 
***************************************************************************************
cd ..
pwd
use "Data\IVR_Data.dta", clear
cd "Output"

drop if endline_completed==0

gen dif_e1_item1 	= e1_cog_1 
gen dif_e1_item2 	= e1_cog_2 
gen dif_e1_item3 	= e1_cog_3 
gen dif_e1_item4 	= e1_cog_4 
gen dif_e1_item5 	= e1_cog_9 
gen dif_e1_item6 	= e1_cog_10 
gen dif_e1_item7 	= e1_cog_11 
gen dif_e1_item8 	= e1_cog_12 
gen dif_e1_item9 	= e1_cog_13 
gen dif_e1_item10 	= e1_cog_14 
gen dif_e1_item11 	= e1_cog_15 
gen dif_e1_item12 	= e1_cog_16 
gen dif_e1_item13 	= e1_cog_17 
gen dif_e1_item14 	= e1_cog_18 
gen dif_e1_item15 	= e1_cog_19 

foreach var in dif_e1* {
recode `var' (0=0) (.=.) (1/10=1)
}

irt 2pl dif_e1_i*, group(any_treat)
estimates store constrained_e1

*item 1
irt (0: 2pl dif_e1_item1) (1: 2pl dif_e1_item1) (2pl dif_e1_item2-dif_e1_item15), group(any_treat)
estimates store item1_e1
irtgraph icc dif_e1_item1, bcc name(i1_e1, replace) title(Question 1) legend(order(2 "Treatment" 1 "Control"))

*item 2
irt (0: 2pl dif_e1_item2) (1: 2pl dif_e1_item2) (2pl dif_e1_item1 dif_e1_item3-dif_e1_item15), group(any_treat)
estimates store item2_e1
irtgraph icc dif_e1_item2, bcc name(i2_e1, replace) title(Question 2) legend(order(2 "Treatment" 1 "Control"))  

*item 3
irt (0: 2pl dif_e1_item3) (1: 2pl dif_e1_item3) (2pl dif_e1_item1-dif_e1_item2 dif_e1_item4-dif_e1_item15), group(any_treat)
estimates store item3_e1
irtgraph icc dif_e1_item3, bcc name(i3_e1, replace) title(Question 3) legend(order(2 "Treatment" 1 "Control"))

*item 4
irt (0: 2pl dif_e1_item4) (1: 2pl dif_e1_item4) (2pl dif_e1_item1-dif_e1_item3 dif_e1_item5-dif_e1_item15), group(any_treat)
estimates store item4_e1
irtgraph icc dif_e1_item4, bcc name(i4_e1, replace) title(Question 4) legend(order(2 "Treatment" 1 "Control"))

*item 5
irt (0: 2pl dif_e1_item5) (1: 2pl dif_e1_item5) (2pl dif_e1_item1-dif_e1_item4 dif_e1_item6-dif_e1_item15), group(any_treat)
estimates store item5_e1
irtgraph icc dif_e1_item5, bcc name(i5_e1, replace) title(Question 5) legend(order(2 "Treatment" 1 "Control"))

*item 6
irt (0: 2pl dif_e1_item6) (1: 2pl dif_e1_item6) (2pl dif_e1_item1-dif_e1_item5 dif_e1_item7-dif_e1_item15), group(any_treat)
estimates store item6_e1
irtgraph icc dif_e1_item6, bcc name(i6_e1, replace) title(Question 6) legend(order(2 "Treatment" 1 "Control"))

*item 7
irt (0: 2pl dif_e1_item7) (1: 2pl dif_e1_item7) (2pl dif_e1_item1-dif_e1_item6 dif_e1_item8-dif_e1_item15), group(any_treat)
estimates store item7_e1
irtgraph icc dif_e1_item7, bcc name(i7_e1, replace) title(Question 7) legend(order(2 "Treatment" 1 "Control"))

*item 8
irt (0: 2pl dif_e1_item8) (1: 2pl dif_e1_item8) (2pl dif_e1_item1-dif_e1_item7 dif_e1_item9-dif_e1_item15), group(any_treat)
estimates store item8_e1
irtgraph icc dif_e1_item8, bcc name(i8_e1, replace) title(Question 8) legend(order(2 "Treatment" 1 "Control"))

*item 9
irt (0: 2pl dif_e1_item9) (1: 2pl dif_e1_item9) (2pl dif_e1_item1-dif_e1_item8 dif_e1_item10-dif_e1_item15), group(any_treat)
estimates store item9_e1
irtgraph icc dif_e1_item9, bcc name(i9_e1, replace) title(Question 9) legend(order(2 "Treatment" 1 "Control"))

*item 10
irt (0: 2pl dif_e1_item10) (1: 2pl dif_e1_item10) (2pl dif_e1_item1-dif_e1_item9 dif_e1_item11-dif_e1_item15), group(any_treat)
estimates store item10_e1
irtgraph icc dif_e1_item10, bcc name(i10_e1, replace) title(Question 10) legend(order(2 "Treatment" 1 "Control"))

*item 11
irt (0: 2pl dif_e1_item11) (1: 2pl dif_e1_item11) (2pl dif_e1_item1-dif_e1_item10 dif_e1_item12-dif_e1_item15), group(any_treat)
estimates store item11_e1
irtgraph icc dif_e1_item11, bcc name(i11_e1, replace) title(Question 11) legend(order(2 "Treatment" 1 "Control"))

*item 12
irt (0: 2pl dif_e1_item12) (1: 2pl dif_e1_item12) (2pl dif_e1_item1-dif_e1_item11 dif_e1_item13-dif_e1_item15), group(any_treat)
estimates store item12_e1
irtgraph icc dif_e1_item12, bcc name(i12_e1, replace) title(Question 12) legend(order(2 "Treatment" 1 "Control"))

*item 13
irt (0: 2pl dif_e1_item13) (1: 2pl dif_e1_item13) (2pl dif_e1_item1-dif_e1_item12 dif_e1_item14-dif_e1_item15), group(any_treat)
estimates store item13_e1
irtgraph icc dif_e1_item13, bcc name(i13_e1, replace) title(Question 13) legend(order(2 "Treatment" 1 "Control"))

*item 14
irt (0: 2pl dif_e1_item14) (1: 2pl dif_e1_item14) (2pl dif_e1_item1-dif_e1_item13 dif_e1_item15), group(any_treat)
estimates store item14_e1 
irtgraph icc dif_e1_item14, bcc name(i14_e1, replace) title(Question 14) legend(order(2 "Treatment" 1 "Control"))

*item 15
irt (0: 2pl dif_e1_item15) (1: 2pl dif_e1_item15) (2pl dif_e1_item1-dif_e1_item14), group(any_treat)
estimates store item15_e1
irtgraph icc dif_e1_item15, bcc name(i15_e1, replace) title(Question 15) legend(order(2 "Treatment" 1 "Control"))

grc1leg2 i1_e1 i2_e1 i3_e1 i4_e1 i5_e1 i6_e1 i7_e1 i8_e1 i9_e1 i10_e1 i11_e1 i12_e1 i13_e1 i14_e1 i15_e1, ycommon labsize(small) position(4) ring(0) lrows(2) cols(4) xsize(8) ysize(8)  scale(.8) name(Figure_B8, replace)
graph export "Figure B8.emf", replace

clear 